/*
  File:    bessel.cpp
  Purpose: Gives any needed Bessel functions. 

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cmath>
using namespace std;

#include "../ImprovedMath.h"

double BesselI( int order, double argument )
{
 double Output = 0.0;
 double CurrentTerm;
 
 int MaxIterations = 20;
 int Iteration = 0;
 double StopAfter = 1e-8;
 
 double ArgSquaredOverFour = square( argument ) * 0.25;
  
 bool stop = false;
 double PreviousTerm = 99;
 
 while( !stop )
 {
  CurrentTerm = IntPow( ArgSquaredOverFour , Iteration ) / 
  ( ( factorial(Iteration) ) * ( GammaFunction(Iteration+order+1) ) );
  Output += CurrentTerm;
  
  PreviousTerm = CurrentTerm;
  
  if( fabs( CurrentTerm ) <= StopAfter )
  { stop = true; }
 
  Iteration++;
  if( Iteration > MaxIterations )
  { stop = true; }
 }

 Output *= IntPow( 0.5*argument , order );
 return Output;
}

